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ABSTRACT 


The break-up of a field of eddies by a flat-plate obstacle embedded 
in a boundary layer is studied using numerical solutions to the two- 
dimensional Navier-Stokes equations. The flow is taken to be 
incompressible and unsteady. 

The flow field is initiated from rest. A train of eddies of 
predetermined size and strength are swept into the computational domain 
upstream of the plate. The undisturbed velocity profile is given by the 
Blasius solution. The disturbance vorticity generated at the plate and 
wall, plus that introduced with the eddies, mix with the background 
vorticity and is transported throughout the entire flow. 

All quantities are scaled by the plate length, the undisturbed 
free— stream velocity, and the fluid kinematic viscosity. The Reynolds 
number is 1000, the Blasius boundary layer thickness is 2.0, and the plate 
is positioned a distance of 1.0 above the wall. The computational domain 
is four units high and sixteen units long. 

A hybrid solution method is used for the velocity field. The 
velocity induction law is used to determine boundary velocities along the 
solid surfaces and on the perimeter of the computational domain. Nonzero 
tangential velocities at solid surfaces are cancelled by the proper amount 

of vorticity production. The velocity field inside the domain is computed 
from the streamf unction. 
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Results are presented over the range of time 0 to 14.61. Vorticity 
contour plots are used to visualize the eddy break-up. Marker particles 
are also used to help visualize the overall flow. A plot of total drag 
variation with time is also given. 

Results show that the eddies are broken up by the plate. The 
strong wake generated by the plate prevents the eddy vorticity from 
penetrating the region between the plate and wall as the eddies are swept 
downstream. Transverse velocities evident ahead of the plate are absent 
behind the plate. Thus, it appears that high speed outer fluid is 
prevented by the plate from being entrained into the fluid layer near the 
wall. This has been proposed to be one of the mechanisms by which break-up 
devices can reduce drag locally. The numerical predictions support this 
proposal. 
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i. introdu ction 


Interest in reducing the drag of aerodynamic surfaces has led to a 
new examination of the fundamental transport processes occuring in 
turbulent boundary layers. It has been hypothesized that the 
proper management of the large scale turbulence can affect wall variables 
such as skin friction. Of particular interest are the large eddy-like 
structures which are believed to entrain the high speed outer potential 
flow into the boundary, thereby causing momentarily high velocities to 
occur near the wall. These higher velocities in turn lead to locally high 
values of the skin friction. 

Whereas there is not as yet universal acceptance of the concept 
that large-scale structures are associated with high skin friction, recent 
important experiments by Nagib and co-workers at the Illinois Institute of 
Technology lend strong support to this concept. Corke, Nagib, and 
Guezennec [1] have found that the outer scales, defined by the intermittent 
excursions of potential fluid into the boundary layer, can be suppressed by 
simple arrangements of parallel plates. This results in a decrease in the 
streamwise growth of the boundary-layer thickness (also the momentum 
thickness), leading to a decrease in the local wall shear stress. They 
report a 30% decrease in local skin-friction coefficient. When account is 
made of the viscous drag of the plates, the overall drag reduction is 20%. 

The authors of Ref. 1 speculate that the large outer scales are 
remnants of the laminar-turbulent transition process. These are slowly 
decaying structures that are embedded in the boundary layer. The plates 
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1 Schematic Diagram of Flow Configuration 



mechnically suppress potential fluid entrainment and thus hasten the 
Reynolds number aging of the boundary layer. The result is a reduced 
growth rate of the boundary-layer thickness with downstream distance. 

The essentially two-dimensional dynamical interpretation given by 
Corke et al. may be an oversimplification of a much more complicated three- 
dimensional flow. However, the basic concept appears to be sound. Four 
mechanisms, by which the plates act to suppress the outer scales, have been 
given by the authors. These can be summarized as follows: 1) restriction 
of the vertical velocity components in the boundary layer, 2) generation 
by the plate of unsteady circulation opposite to that of the large-scale 
motions, 3) generation of a small-scale vortex street behind the plate, 
and a redistribution of small-scale turbulence, and 4) small-scale 
turbulence production due to the wake of the plates embedded in a much 
thicker wall boundary layer. These above-mentioned mechanisms are not all 
independent, and when account is taken of the unsteady nature of the flow 
and the many scales of turbulence present, it is difficult to identify 
which ones are dominant. 

The present study was undertaken in order to identify the basic 
dynamical processes which occur when eddy— like structures interact with a 
flat plate embedded in a boundary layer. The approach is numerical rather 
than experimental and is based on solutions to the unsteady two-dimensional 
Navier-Stokes equations for incompressible flow. The intent is not to 
simulate a fully turbulent flow by numerical means, nor is any turbulence 
modeling incorporated into the analysis. Rather, an analysis is made of a 
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flow which has several gi^oss features in common with the boundary-layer 

flow studied experimentally in [1]. 

The eddy-like structures are taken to be regions of constant 
vorticity which are initially circular in shape* They distort in the shear 
flow as they are convected towards a single flat-plate manipulator 
positioned parallel to a wall of large extent* A sequence or train of 
eddies is introduced computationally ahead of the plate. The eddies have 
varying strengths and length scales, the latter being on the order of the 
boundary- layer thickness, <S . These length scales are comparable to those 
found in the experiments of [1]. 

The eddies are superimposed onto an otherwise steady background 
flow which is laminar. In this way, a flow is produced which is close to 
that in a laminar boundary layer undergoing transition. Thus, it is not as 
developed as that produced experimentally in [1]. Also, the Reynolds 
number based on the boundary- lay er thickness is 2000, which is 
approximately 1600 when based on the momentum boundary-layer thickness. 
This latter quantity ranged from about 2200 to 5000 in the experiments of 
[l]. Therefore, this parameter for the present study is close to the low 
end of the range studied in the experiments. 
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II. BASIS OF THEORETICAL APPROACH 


The basic approach is the same as that incorporated by Schmall and 
Kinney [2]. However, rather than the flat plate being positioned in a 
uniform onset flow, in this work it is imbedded in an unsteady viscous flow 
adjacent to a plane surface. This necessitates the introduction of image 
vorticity in order to satisfy boundary conditions. As in [2], the plate is 
replaced by a distribution of bound vorticity, and the vorticity production 
at solid surfaces (plate plus lower bounding wall) is calculated directly. 

The velocity field is calculated using a hybrid scheme which 
incorporates the stream function and velocity induction law. The velocity 
induction law is used to obtain the effect of the bound vorticity of the 
plate on the main flow, as well as the tangential velocity induced at solid 
surfaces. This latter quantity is needed in the determination of the 
vorticity production. The rotational velocity of the main flow is computed 
from the stream function, as obtained from the numerical solution of 
Poisson's equation. 

The analytical development is given in the next section. The 
working equations are presented, and the numerical formulation is 
discussed. Lengthy developments are given in the various appendices. 
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III. ANALYSIS 


A schematic of the flow configuration is shown in Fig. 1. The 
characteristic length scale for the flow is taken to be the plate 
length, L. We have selected S = L and 5 = 2L. The background velocity 

P r °fil6 is denoted by U^(y). This is taken to be the Blasius profile for 
laminar flow on a flat plate with zero pressure gradient. 

Strictly speaking, 6 and are functions of x. However, 

over the region of interest, their streamwise variation is very small. 
This can be seen if the characteristic Reynolds number is introduced as 
“ UqL/v where a value of 1000 has been chosen for this study. 

Let x 1 be the streamwise distance measured from the virtual origin 
of the boundary layer on the wall. One has for a Blasius profile that 
6/x' = b/CU^xYv) ^ . Then <5/L = 6[ (x'/W/CU^L/v . Upon substitution 
of 5/L = 2 and U^L/v = 1000, one obtains x*/L = 111. Therefore, the 
virtual origin of the boundary layer is more than 100 plate-lengths 
upstream of the plate. As will be discussed more fully later, the domain 
of interest spans a region which is 16 plate— lengths in the flow direction. 
Over this distance, the boundary layer grows only 7%. This is a low order 

effect, and therefore the variation of 5 and U^ in the streamwise 
direction has been neglected entirely. 

The other quantities of interest are the vorticity of the 
background flow and the relative distance y/L, where y is 
measured from the wall. One begins with the boundary— layer approxi- 
mation 9y, where now is considered to be a function of n. , 
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and n = (y /x')(U x'/v) Recall that 0 = U f'(r|), where f'(tl) 

oo D 00 

is the derivative of the dimensionless stream function. Clearly, 

- 1/2 

fl --U^f'Cn) 3n/8y, from which O b L/U 00 = -f'Cn) [(x'/L)/(U C0 L/v)] is the 

dimensionless background vorticity. From the foregoing expression for 

<5/L, one has [(x'/L)/(U KL/v)]^ 2 = 1/3. Therefore, ^ b L/U OT = -3 f"(n). 

- 1/2 

One can also show that n = (y /L) [x'/L) (U^ L/v)] , from which 

y/L = n/3. When H= 6, then y = 6 and <5/L = 2 as required. 

To summarize, then, we have 


°b - 

V*(n) 

(l) 

« b L/u » * 

-3 f"(n) 

(2) 

y/L 

n/3 

(3) 


The values of 0, f'Ol), and f"(il) are given in tabular form in Ref [ 3 1 » 
Page 139. In practice, it is the value of y/L which is specified in the 
calculations. From this, q is obtained from (3), and the values of f'Ol) 
and f"(Ti) are obtained by interpolation. 

Governing Equations 

Superimposed onto the background flow are perturbation velocity and 
vorticity fields caused by the presence of the plate, as well as the 
coherent eddy structures introduced ahead of the plate. We denote as u' 
and v' the x- and y-components of perturbation velocity, and or' denotes the 
perturbation vorticity field. Then 
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( 4 ) 


u = U b^ y ^ + 

v = v' (x,y,t) (5) 

0) = ^ b (y) + w'(x,y,t) (6) 

No assumption is made that these perturbation quantities are small. 

The vorticity transport equation is given by 


3co 

3t 


+ 


s (u “ J + 



_1 (IJ* 

Re 3y 2 


3jo N 

3x 2y 


(7) 


and the continuity equation is 


3u 9v 

3x 3y 


0 


(8) 


At this point, all independent and dependent variables have been rendered 
nondimensional using U^, L, and v. That is, t ~ t^D^/L, u = u^/U^, y = 
y*/L, to = to^L/U^, etc. The asterisk denotes a dimensional variable, which 
will normally be omitted for convenience. The origin of the x and y 
coordinates is as shown in Fig. 1. 

Note that Eqn. (7) allows the background vorticity, , to be 
transported by convection and diffusion in the lateral direction only, 
since it does not depend on x. Therefore, there is a mechanism by which 
can change with time. Nevertheless, the effect is small and is of no real 
importance. Therefore, it has been neglected here. 

The velocity field is obtained from a hybrid formulation in terms 
of the stream function and/or the velocity induction law, referred to here 
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as the Biot-Savart law. The general form for the x-component can be 
written as 

! ff “ Cy-y ) 

uCx p’y p t) = 2? , 7z , 72 dx dy + i*g rad ^ (9) 

11 (x-x ) + Cy-y ) 

where P is some point in the flow, and the range of integration is over 
the region of non-zero vorticity. Note that this expression contains 
i ‘grade}), which must be included for generality. More will be said about 
this term later. The expression for the y-component of velocity is 

1 (■ (• w (x-x ) 

v (x , y ,t) = - - 77 - r 2 dx dy + j ‘grade}) (10) 

p p JJ Cx-x p ) 2 + Cy-y p ) 

It is easily verified that curl v = k oj , where the curl operator is with 

respect to the coordinate at P. 

The term grad $ is a purely irrotational contribution to the 
velocity field and must be included to insure that boundary conditions are 
satisfied. The principal velocity boundary conditions embody the adherence 
condition at solid surfaces. These are enforced in two steps. First, the 
normal velocity component is nullified, after which the tangential 
component is reduced to zero. The first step is accomplished through image 
vorticity plus the proper specification of grad<j>. The second step is 
accomplished through the proper production of vorticity at the solid 
surfaces . 

We require that v vanish on the solid wall and plate. The wall is 
taken to be a plane of anti-symmetry such that below it we have vorticity 
which is opposite in sign to that above the wall. This concept of image 
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vorticity provides a straightforward means for determining grad <p in (9) 
and (10). 


Any irrotational velocity field which produces zero normal velocity 
at the plate and wall will satisfy the requirements for grad<J>. Therefore, 
we replace the plate by a distribution of bound vorticity, y , along its 
length. A distribution of fluid sources would do as well, but the 
vorticity distribution better satisfies our needs. Therefore, we make the 
following substitution. 


grad <p - — ± ~2~ dx (11) 

- 1/2 r 

where r is the vector from y to the point P, and dx is the incremental 


plate length over which y is distributed. Also, y = k y. 

To produce the proper anti-symmetry about the wall, we replace co 
and y in the lower half plane (y * < 0) by their images. Therefore, we 


have 


U( VV C) = U 


4* OO — J— oo 


■ f y - y o 

J ^ (x-x ) 2 + (y-y ) 2 

r\ — co r> w 


y-y D , y+y P 1 , . 

2 + 2 2 dx dy 

+ (y-y ) 2 (x-x ) 2 + (y+y ) Z J 


t 1/2 r -7 

Y J £ — — J 

i ,/, + ( s -yj 


2 2 dx 

(x-x p ) + (s+y p ) J 


OO 4* OO 


v (x »y n ,t) = 
p p 


0 — oo 


-x p ) 2 + (y-y p ) 2 


J dX dy 

(x-Xp) + (y+y p ) J 


+1/2 


-x p ) + (s-y p ) ( x_x p ^~ + (s+y p )' 



To obtain (12) and (13), one replaces co with ~co and Y with -y in 

the lower half plane. The range of integration on y is then changed to 

include only the positive half-plane y >0. 

It is clear from (13) that v = 0 whenever y =0. That is, the 

P 

wall is a streamline. Note, however, that u ^ 0 when y p = 0 , and there 
is an apparent slip velocity on the wall (there is also one on the plate). 
This must be reduced to zero by the proper production of free vorticity at 
the wall and plate. This process will be discussed subsequently. 

The distribution for Y must be so constituted that no transverse 
velocity exists at the plate. To obtain the governing equation for y , one 
sets y p = S in (13), and sets the left-hand side to zero. There results 
an integral equation for y as follows: 

(14) 


+1/2 +03 J 

h 00 


f -1 X-X 

yf 1 P 1 dx _ 


, 

x-x X-X -I 

P P 

Y x-x . .2 2 dx 

-1/2 L p (x "V + 4s J ( 

J 

) - 

- 00 

_(X-Xp) 2 + Cy-y p ) 2 (x-x p ) 2 Hy+yJ 


When it is solved, the resulting velocity field satisfies all the 
requirements of the problem, except possibly the no-slip condition on the 
wall and plate. As mentioned earlier, this must be reduced to zero in the 
calculation of the free vorticity from (7). 

For the evaluation of the apparent slip velocities at solid 
surfaces, one uses (12). It is expedient to introduce co = ft^(y) + 
<u'(x,y,t) into (12) and to note that the integration of ^(y) over the 
infinite domain produces only the background velocity, ^b(y p )* Therefore, 
we add N^(yp) to the right-hand side of (12) and replace co by co *. Next, 
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t?e set y = 0 and evaluate u(x ,0,t) . Let us call this (U - . ) 


slip wall 


Thus , 

slip' wall 


1 f" f" 2y 

- 2tT j j T 1 dxdy 

0 


(x-x p ) + (y-y p )‘ 


(15) 


1/2 




2S -j 

2tt 

Y 

-1/2 

(x-x p ) + (S-y p ) 2 


dx 


The evaluation of the slip velocity on the plate is slightly more involved. 
First, set yp -S + in (12), where S is a small distance £ above the 
plate. The integration of co is straightforward. The integration of 
the first term involving Y must be done carefully, since y = S , and thus 
S - y p = -e. As e + 0 we obtain -j/2 for this integral. When 
y p a S “ we obtain +y/ 2 for this integral. Thus there is a different 
apparent slip velocity on the upper and lower surfaces of the plate. The 
end result is 




^slip^plate 


U b^ + 2^ 


-f OO 4-00 

f f 


0 - 00 


(o' ^(x.y.Xp.Yp) dx dy 


(16) 


+1/2 


1 , N . 1 

- 2 ^ ( V + 2i 


Y(x) 


- 1/2 


- 2S 

2 2 
(x-x V + 4S 


dx 


(D slip ) plate ’ (B .lip>plate + y( V 


( 17 ) 


where K is the quantity in brackets in (12). 
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Once the slip velocities are determined for the wall and plate 


surfaces, the vorticity production is given by 
t+ At 

f — ( dt* = (TJ ) 

J Re \9y /surface v slip surface 

t 


(18) 


This follows the scheme used by Schmall and Kinney [2], The convention is 
to evaluate the right-hand side of (18) at time t. Also, a positive slip 
velocity on the upper portion of a surface produces negative free 
vorticity, whereas the reverse is true on the lower portion of a surface. 

Equation (18) is the essential surface boundary condition needed in 
the solution of (7). The manner in which it is implemented is discussed in 
the next section. 

As mentioned earlier, a hybrid formulation was finally used for the 
calculation of the velocity field at interior points of the flow. This 
involved the use of the stream function for a part of the evaluations, as 
is now explained. 

The integral involving co, or equivalently o', in (12) can be cast 
in terms of the stream function. All that is required is that this 
integral give a rotational velocity field which has the correct amount of 
vorticity at each fluid point and which satisfies the boundary conditions. 
That is, we can replace this integral with U^(y^) + / 3y , where 4 r ' i s 
the perturbation stream function which satisfies the Poisson equation 


2 ' 2 

w + o: 
2 2 
9x Sy 


-co 


(19) 
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It turns out that the evaluation of dfy'/dy is computationally more 
efficient than the evaluation of the integral involving oj ' in (12), when 
the domain of integration is large. It remains to ensure that the proper 
boundary conditions are specified for ifj' or its derivatives. 

We first write 


+ CO -f 00 


= _1_ 

3y 2 tt 


1 

3x 2ir 


co' Kj-Cx.y.Xp.yp) dx dy 


0 - 00 
4-oo 4” 00 


0 - 00 


«' K TT (x,y,x ,y ) dx dy 
11 P P 


( 20 ) 


( 21 ) 


where K^. and K^. are the quantitites in brackets in (12) and (13), 
respectively . 

For a given w 1 , (20) or (21) give the value of the normal 

derivatives of ijr* at points on the boundary of the domain, inside of which 
\fj* is defined by (19). That is, if y^ is set to the value 4.0, then (20) 
gives the value of dipVdy, at points x^ along the top (horizontal) boundary 
of the domain. Similarly, for fixed x^, (21) gives the value of 3^*/3x at 
points y^ along the side (vertical) boundaries of the domain. On the wall, 
we set =0, since this must be a streamline. 

The boundary condition along the wall is of the Dirichlet type, 
whereas those along the sides and top are of the Neumann type. It is known 
that the solution to (19), subject to these boundary conditions, is unique. 

Before turning to the numerical formulation, a few comments 
concerning the evaluation of y are i n order . Recall that the integral 
equation for y is given by (14). The vorticity field, 0) 1 , is known. The 
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solution for y induces a purely potential (i.e., irrotational) velocity 
field, which corresponds to that produced by a flat plate in ground effect. 

this case, the velocity induced at the plate by co* is simply a given 
input quantity equivalent to an outer onset-flow velocity. 

As S the solution for y must approach that for a single 

plate in unconfined flow. It is equivalent to that given in classical 
aerodynamics and used by Schmall and Kinney [2]. Therefore, the general 
form is already known from classical theory. 

The essential point to keep in mind here is that the solution for 
Y is non— unique and must contain a term corresponding to pure circulatory 
flow about the plate. For a single plate, this term is A/sin 8. Here A 
is a constant and 0 is a polar angle obtained by circumscribing a circle 
about the plate so that the plate is coincident with a diameter. The angle 
6 is measured from the plate in the counterclockwise direction. For a 
plate in ground effect (i.e., there are two plates when the image is 
introduced), the solution corresponding to the pure circulatory flow is of 
the form A f(9)/sin0. Clearly, f(0) -► 1 as S + °°. Unfortunately, 
f(0) is not given in terms of elementary functions. However, it is easily 
determined numerically. It is found by solving the integral equation for 
y when the non-homogeneous term is omitted. The non-homo geneous term 
corresponds to the velocity induced by cu ’ . The resulting solution for y 
to the homogeneous equation is called the complementary solution. 

The unknown constant A is determined from the principle of 
conservation of total vorticity (see Ref. [2], Simply stated. 
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( 22 ) 


+ 1/2 ^ /N 

/ Y dx = j Y C.9)d6 = 0 

- 1/2 0 

/s ~ 

where Y = Y /sin0. Recalling that Y (or Y ) is composed of a particular 

solution to the complete (i.e. non-homogeneous) integral equation plus the 

complementary solution [i.e. Y = y + A f(0)], A can be determined 

pare 


to be 


TT 


- / Y „ (0) de 

1 part 


!*" 

0 


f C0) dQ (23) 

where it is convenient to carry out the integrations in the 0-plane. 


Numerical Formulation 

A complete description of the computational domain is given in the 
section dealing with the numerical parameters. Briefly, it is rectangular 
in shape. It extends from x = -4 to x = +12, and y = 0 to y = 4. A view 

of the region near the plate (-4 - x - 4) is given in Fig. 2. The plate is 

< < 

positioned at y = 1, and it covers the range -1/2 - x - +1/2. The 
background velocity field is confined to the region 0 - y - 2. That is, 5 
in Fig. 1 is equal to 2.0. 

It is assumed that ca 1 is zero above the computational domain y > 4. 
Furthermore, a) 1 is set equal to zero ahead of the upstream boundary 
(x < -4) and behind the downstream boundary (x > 12). At x - 12 we have 
imposed the condition 9a) f /3x = 0. 

The grid is refined near the wall and plate in the transverse 
direction, and near the leading and trailing edges of the plate in the 
streamwise direction. Furthermore, the plate is discretized along its 
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O 
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DISTANCE ALONG THE WALL, X 


Fig ♦ 2 Computational Domain Near the Plate 


length independently of the outer— flow grid. There are 20 fluid cells in 
contact with the plate, whereas there are 80 points along the plate at 
which y is determined. 

Solution for Bound Vorticitv. The solution for y is obtained at 
80 discrete points on the surface of the plate. The polar angle 0 is 
introduced, and the 80 points are obtained by first constructing arcs of 
length A6 on the circle of unit diameter, where A0 = tt / 80. Points are 
positioned at the centers of each arc around the circle, and these are 
projected vertically downward onto the diameter, which comprises the plate. 
In this way, there is a clustering of vortex points on the plate near the 
leading and trailing edges. 

A polar coordinate system is used in which x *= .5 cos0 and 

dx ® - .5 sin0d0. The integral equation for y, Eqn. (14), is evaluated 

at each of 80 points given by x = .5 cos 0 • The solution is actually 

P P 

^ A 

determined for y ( = ysin0) rather than y , and y is assumed to be 
constant over an interval A0. When this is substituted into each of the 

80 integral equations, there results a system of 80 simultaneous equations 

/% 

in 80 unknowns, y^. 

As previously mentioned, the general solution is composed of 
particular and complementary solutions. Since each is non-unique, we can 
assign an arbitrary value to one of the 80 unknowns in each solution. 

For convenience, we set the 80th unknown in each set to 1.0 and 
solve the two sets of equations. If f(0) denotes the solution to the 

A 

homogeneous equations, then Y com p = A f(0). The constant A is 
determined from Eqn. (23). 
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Solution for Free Vorticity « The solution of Eqn. (7) is 
straightforward. We first integrate (7) with respect to space over a fluid 
element AxAy. This gives 


3t 


dx dy = - 


+ te 


div [i u to + j v co] dx dy 
div (grad gj) dx dy 


(24) 


By the divergence theorem. 


y+Ay 


div [i u u + j v 4)J dx dy = 
x+Ax 

+ 


[C-uco) x + (U£0) x+Ax ] dy 


(25) 


[C-vu) y + Cvw) y+Ay ] dx 


and 


div (grad w) dx dy = 


y+Ay 

[ 


y 

x+Ax 


/ 3oo\ , / 3co 
\~ 9x/x \3x 


x+Ax 


] dy 


(26) 




We apply the mean-value theorem, which leads to the approximation 

y+Ay 

(£)... * ! (g) ^ 


x+Ax 


( 27 ) 


x+Ax 
y+Ay/ 2 

and so forth for each of the terms in (25) and (26). If subscripts i and 
j denote the variables evaluated at points in the vicinity of a fluid 
cell, as shown in Fig. 3, then we can write further that 


y+Ay 


/ jko \ , ~ r ^i+l.i - ^i.j ~l . 

( 3x ) *+Ax y _ [ ' 5 C4x i +4x i+l> ] A 


( 28 ) 


3t 


£0 dx dy = -r— [w. . ] Ax. Ay. 

9t i,j i J 2 


(29) 
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Fig. 3 Diagram of Node Arrangement around Fluid Cells 
in the Computational Domain 


I I 



and so on and so forth. Note that the nodes are midway between control 
faces (dashed lines). For the evaluation of the convective terms in (25), 
we use an upstream weighted scheme as follows: 


y+ A y 

(uco) 


x+Ax d ? ' t“i+l/2, 3 “i.j 1 % 


(30) 


Note that the velocity at the right-hand cell face convects vorticity at 
the upwind node. 

When these expressions are assembled into (24) and each term is 
divided by Ax_^Ay_. , the result is 

“i,j = £r{ [u i-l/2,j Vl,j ] ■ [u i+l/2,j 


Ayj { [V i,j-l/2 “i.j-l 1 [V ±,j+l/2 “i.j 1 ) 

2 < _1_ r (0J i+l ? j ~ 

j Ax ± LCAx i + Ax i+1 ) (Ax ± _i + Ax ± ) J 


Re 


r (u i,j + i - 

CO. .) 

i,.l 

(co . . — (0 . . . ) -] X 

1,3 i,3-l 1 


L(i Vi + 

Ay t > 

(Ay., + Ay. )Jj 

(31) 

i integrate 

(31) 

with respect to time. 

We use the 


explicit method which presumes that the right-hand side of (31) is 
evaluated at time level t . Next, both sides of (31) are integrated over 
time interval At . Upon rearrangement, one has 
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1 


CO 


k+1 = w k . . / i-u k AL _ k At _ 2 At r 1 / 1 , 

*•>3 ( i+l/2,j Ax. i,j+l/2 Ay. Re Ax. lAx.+Ax.,- Ax 4 , t +Ax. 

v i 3 L x \ i x+1 i— 1 x 

-1 — / i + i Yh+/,\k r 2At/Re — 

4y j-i +4y j)|} i+i>j L Ax i (4 v Ax i + i>] 

T 2At/Re 1 

h wy j + 4y i + i ,,j "j 


JL ^ 

+ Vi.j 

+ ftf 
l.J-1 


r k 


Ax. 

i 


2/Re 

Ax. (Ax. .+ Ax. ) 
x x-l 1 


At + <0 


i»j+l 


Ui-i/2 

A __ * 


2/Re 


At 


Ay.CAy._ 1+ Ay.)J (32) 

Note that the superscript k now denotes the time level. 

For nodes adjacent to solid surfaces, Eqn. (32) must be modified to 

account for zero normal velocity and the vorticity production. One sets 

v. . , = 0 or v. = 0 , depending on whether the cell is bounded 

x,j-l/2 x»j+l/2 

below or above by a solid surface. The vorticity production term replaces 
the diffusion term given by 


t+At x+Ax 


\ \ -(t). 


dxj dt 


(33) 


x 


where either y = 0, y = S , or y = S . This entire term is replaced 


v+ 


by - (D , . ) Ax in the case of y = 0, and ~ (U , . )' Ax or 

slxp wall J slip plate 


(U , . ) , Ax in the case of y = S 

slxp plate 


or S , respectively. 
Accordingly, one sees that (32) must be modified by adding 


(co^ . - a) k . n ) 2 At/Re Ay. (Ay. , + Ay.) 
x,j x,j-l 3-1 J 

to the right-hand side. This zeros out the flux terms at the lower surface 
of the control volume. Finally, - (E^^ ) sur f ace /Ay > is added to the right- 
hand side [recall that (33) must be divided by Ax.Ay^]. A corresponding 
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modification is made for y = S , 
added. 


except that now +(tJ , . ) , „ /Ay, 
r slip plate J 1 


is 


Observe that initially k = 0 , and the vorticity is everywhere 
zero. One then has the simple result that adjacent to the top of solid 
boundaries , 




slip 


ip 


Ay, 


(34) 


which shows that the initial slip velocity, which after all corresponds to 

a vortex sheet on the surface, produces the correct amount of free 

vorticity at the node i,j in time At. That is, the initial sheet vortex 

of strength (U ) is broadened by diffusion in time At until it fills 
slip i 

the whole fluid cell adjacent to the solid surface. Therefore, At must be 
chosen to be the diffusion time for a fluid cell of height Ay. This also 
ensures the stability of (32), as explained below. 

Equation (32) will be stable for sufficiently small At such that 
the coefficient of the a). . term is positive or zero. Adjacent to the 
solid boundaries, the convective velocities can be neglected, as can the 
streamwise diffusion of vorticity. Then (32) will be stable provided 


1 


2At 

ReAy^ 


Ay. + Ay 


— ) > o 
j+VJ 


(35) 


This is precisely the statement that At be on the order of the time for 
for vorticity to diffuse a distance Ay. Tar from solid surfaces, the 
incremental heights Ay are greater, and the fluid velocities are 
important. Then At must be sufficiently small that the entire coefficient 
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k • 

of the co term in (32) be positive or zero- In all cases, the stability 
■L* 3 

critereon was satisfied- It was also verified that the Ay spacing next to 
solid surfaces was on the order of the diffusion distance for time At* 

Solution for Stream Function and Velocity Field, As will be 
explained more fully later, the solution for the perturbation free 
vorticity, co 1 , is obtained first. Therefore, the vorticity field is always 
considered to be known at any given time level. 

The solution for the stream function follows from Eqn. (19). For 
the nodal arrangement shown in Fig. 3, one has by the circulation theorem 


$4 , +1 ~ % A 

- co. . Ax. Ay. = - —t * 3 . Ax. 

i J .5 (Ay. + Ay ) i 


- *± ,1 - Ay . *±,i - ♦i. j-i Ax 

.5(Ax ± + Ax ±+1 - y j •5(Ay_._ 1 + Ay^) i 

ty. , . i|>. . 

4 1 A? j Z ill A v 

+ . 5 (Ax . . + Ax.) Ay j 

i—l i J 

where the primes on 1 p and to have been omitted for convenience. 

Upon rearrangement, an equation for ib can be written as follows: 

-1 


(36) 




itj-1 


where B = 


B co. . + C . + D ill. .,,+E xh . , _ . + F , . 
i»J i.J v i,j+l Y i+l,J v i-l,j 

Ay. 

-i- + a v 


(37) 

(38) 


Ay. (Ay i+Ay.) (Ay. i+Ay.) Ay (Ay ,+Ay ) 

c = 1 + — 3 lZ± 1 + Hi 1 + — 2 Hi J_ 

Ax.(Ax.+Ax . .) (Ay.+Ay.,.) Ax. (Ax. .,+Ax.) 

i i l+l J 2 j+l i i-l j 

D = - ( Ay i-l +fly i ) 

(A yj +Ay J+1 ) 


(3 9) 
(40) 
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Ay i (Ay^+Ay..) 

Ax. (Ax.+Ax. . n ) 
1 1 i+l 


(41) 


E = 


Ay (Ay ,+Ay.) (42) 

T7 = _ J U — 1 

Ax. (Ax. ,+Ax.) 
l l-l l 

The foregoing equation applies at any interior node for which the control 
volume (dashed lines in Fig. 3) does not lie adjacent to a boundary. A 
special form must be used for boundary nodes, since the boundary conditions 
corresponding to Eqns. (20) and (21) and = 0 on the wall must be 
satisfied. The task is to find values corresponding to ty. , . and 
at left- and right-hand boundaries, respectively, as well as values 
corresponding to above the top boundary. The boundary condition at 

the wall is of a different type and will be discussed subsequently. 

Recall that the boundary conditions on the sides and top of the 
domain are given by Eqn. (20) and (21), where the right-hand sides are 
considered known. Consider first the left-hand boundary. Let (Xp,7p) be 
the node (i-l,j) as shown in Fig. 3. Then 



(43) 


and the quantity v. is known. Now expand ip in the Taylor series as 

L 

follows : 


ip. . 

itJ 


= Tp 


i-1. j 


♦(£) 


Ax + 


1-1*3 


/lijh 


Ax 

2 


0(Ax 3 ) 


(44) 




i+l. j 



i~l* j 



(45) 
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Now eliminate (S ip/dx ). , . from the above to get 

1 -L| J 


®u 


4\p. ... -^ 1 .- 3 ^.,. „ 

— ^ + 0 (Ax2) 


( 46 ) 


Note that we have made use of the fact that adjacent to vertical boundaries 

we have uniform spacing such that Ax. . = Ax. = Ax. ,, = Ax. Now substitute 

l—l i l+l 

(43) into (46) to get 




2Ax v L + - 4» 1+1-1 


(47) 


The corresponding expression on the right-hand boundary is 


- 2Ax v + 4ip - ip 

iji = £ — ixJ A.A? J 

V i+l.j 3 


(48) 


For the top of the boundary we use the simple first-order scheme. 
This is necessary due to the solution algorithm adopted to find . 
This will become apparent shortly. One obtains for the top row of nodes 


’"i.j+i ' '"i.j + u i 4y 


(49) 


where u T is known from the right-hand side of (20) for (x ,yp) 

corresponding to the node (i, j+1) . 

The enforcement of the boundary condition at the wall is treated 

next. It is recognized that the wall coincides with the horizontal 

control-volume face between nodes i, j-1 and i,j. Along this face, ip. 

i,w 

must be zero. We first obtain an expression f or using the procedure 

i,w 
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adopted to derive Eqn. (34). That is, we apply the circulation theorem to 

a new control volume surrounding node i,l. The height of this control 

volume is .75 Ay^ . The node i,l is thus a distance of .5Ay^, above the 

wall and .25 Ay^ above the lower control face. The upper horizontal 

control face is a distance of .5 Ay^ above node i,l. The resulting 

equation differs from (36) in only minor respects. When the result is 

rearranged to give \b one obtains 

i.w 


*i,w = -(• 375i J'l 2 ) 


CO 




1 + 


A yi 


Ay l +Ay 2 


+ . 75Ay ' 


Ax.(Ax.+Ax ) 

• i i i+i 


+ Ax i (Ax i _ 1 +Ax i ) ] { ^i,l ~ (Ay 1 +Ay 2 ) ^i,2 

' 75 Ay i 2 i - 15 Ay i 2 i 

[ax^Ax^+Ax^j. J ^i+1,1 " fix i (Ax i _ 1 + Ax i ) J 'i-1, 1 (50) 


Modifications corresponding to Eqns. (47) and (48) are introduced at the 

upstream and downstream boundaries. As discussed below, it is required 

that be zero for all i. 

r ,w 

The solution algorithm used to calculate ijj, . is based on the 

i > J 

Stabilized Error Vector Propagation (SEVP) method described in Ref. [4 ]. 
The technique is due to Madala [ 5 ] and is an extension of the scheme 
described by Roache [ 6 ] . 

Basically, one begins at the top of the computational domain with 

j=N. Eqn. (49) is used to replace \jj m in (37), and expressions (47) and 

x , j • 1 

(48) are used adjacent to boundaries. The result is an expression for 
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d; „ „ in terms of \ b. „ „ 9 ty. - T ,11;.,- XT and other known quantities. If 
y i,N-l r i-l,N i,N i+l, N 

values for the stream function along this top row of nodes are given 

arbitrarily, then one can extend the solution to the next lower row of 

nodes. This is done until the nodes are reached at the wall, and the 

stream function there is calculated by (50)* Had the initial distribution 

of \b along the top row of nodes been correctly given, then given by 

1 ,w 

(50) would be zero. The departure of from zero is a measure of the 

1|W 

error, and one can correct the guess for ^ along the top row of nodes 
according to a systematic procedure. In this way, the correct distribution 
at j=N can be deduced, and the true values for the stream function over the 
whole field can be calculated. This is a direct solution method which 
avoids iteration. 

As pointed out in [6], the method fails for large computational 
grids. That is, any error introduced at some row will eventually grow 
until the solution becomes meaningless. Therefore, the solution has to be 

stabilized by subdividing the vertical extent of the domain into subregions 

or blocks. A direct measure of the error can be obtained, and so 

subdivisions proceed until the error is acceptable. 

The numerical algorithm actually employed is quite complex and is 

too lengthy to describe here. Details are given in Appendix A. 

Having obtained the stream function, it remains to obtain the 

velocity component, u = 3ip/dy. Since \p is computed at grid points, the 

usual expression given by u = (ip. ... - ip. .)/-5(Ay.+ Ay.,.) gives u 

j-jJ J 

which is tangential to the top of a control volume (see Fig. 3). The value 
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of u desired is perpendicular to vertical control faces between node 
points. Therefore, a different scheme was devised. 

A biquadratic surface was used to interpolate values of ^ between 
nine node points. The polynomial is of the form 


2 2 2 o o o 

^ = a + b x +cy + dxy + ex +fy +gxy + hxy + ix y (51) 


from which 


u = 3if>/3y = c + dx + 2fy + gx + 2hxy + 2ix^ 


(52) 


The nine unknown coefficients (a, b, i) are found in terms of the nine 

values of t|) at node points (i-1, j — 1) , (i,j-l), (i+1, j-1), ..., (i+1, 

j+1). Then u is found from (52). For example, the value of u. , . at 

i-h,j 

the vertical control face between nodes (i-l,j) and (i,j) is found from 
(52) with y = .5 (Ay.^ + Ay ) and x = .SAx.^. Unfortunately, the 
coefficients c, d, f, etc. are not given by simple expressions when the 
node spacing is variable, and therefore, they are not given here. However, 
they are summarized in Appendix B. 

Evaluation of Velocity Field due to Bound and Free Vorticity. The 
presence of the plate in the flow field means that the domain is not simply 
connected, and the flow can have a purely irrotational component. As we 
have seen, this can be handled by the introduction of bound vorticity on 
the plate. One must calculate the velocity field due to this bound 
vorticity. The complete expression for u is given by 
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1/2 

U <V V 0 ■ " b<y) + If + n 1 Y 


s-y. 


s+y. 


- 1/2 


F 


x ) 2 +(S-y ) 2 
P P 


(x-x p ) 2 +(S+y p ) 2 


Idx 


(53) 


The integral over the plate in the above expression is done numer- 
ically. First, polar coordinates are introduced, as discussed earlier 
(i.e. x - .5 cos 8). The quantities x p , y p and S( = l) are fixed. Next 

Y sm0 is replaced by y Over any one of the 80 intervals A0 * v is 

i 1 i 

assumed to be constant. For points close to a plate segment, the quantity 
in brackets in (53) is integrated numerically over the segment using the 
Gauss quadrature method with 20 Gauss points. For points far removed from 
the segment, the quantity in brackets is approximated as being constant 
over a segment. It is evaluated at the midpoint of the segment, and the 
integral is this value multiplied by A9_^. Following this procedure for 
any fluid point (x^, y p ), one obtains a sum of 80 contributions for the 80 
values of . Each contribution is a value of Y^ multiplied by a 
geometrical coefficient, which is constant for all time. Therefore, these 
80 coefficients for each fluid point were calculated once and stored. 

Integrals over the field of free vorticity exist in Eqns. (14), 
(15), (16), (17), and (20) and (21) for (x p ,y p ) on the right, left, and top 
boundaries of the computational domain. Recall that these latter 
expressions are needed to find v R , v^ and for use in (47), (48), and 
(49). In all cases the basic idea is the same. Over any fluid cell, co . . 
is considered to be constant and is factored outside the integral. The 
kernel functions are then integrated over cells Ax^ Ay. , Since each cell 
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is rectangular, it is divided into 40 vertical strips, which are centered 


on the 40 Gauss points for the interval A-x^. Each vertical strip is 
treated as being of infinitesimal width, and the integral over Ay. is 
carried out exactly (see Ref. [ 7 ] for more details). These results are 
then multiplied by the corresponding dx_^, and the results are 
summed according to the weighting scheme given in the Gauss quadrature 
formula. This gives one geometrical coefficient for the cell Ax^Ay . 
and a point (Xp,yp). When the cell is far from the point, that is 
r > 10 ( Ax., Ay ) ' , then the kernel is treated as constant. The integral 

1 j 

is then the kernel evaluated at the center of the cell multiplied by the 
cell area. 

In this way, geometrical coefficients are computed for each fluid 
cell and any given point (Xp,y^). The results are calculated once and 
stored for each point. 

In principle, the y-component of velocity could be found from an 
equation similar to Eqn. (53). It is more expedient to obtain it from the 
continuity equation, once the u velocity component is found. That is, 
Eqn. ( 8 ) is integrated with respect to y and then with respect to x. Thus 


f y 

u(x+Ax,y)dy 


f y 

u(x,y)dy + 


x+Ax 

r 


v(x,y)dx = 0 


o ox 

where the wall condition v( 0)=0 has been applied, 
integrals can be approximated as follows: 


(54) 

To second order, the 
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0 


(55) 


n n a 

E u (x+Ax , , y . ) Ay . - E u(x,y.)Ay. + v(x+-r— , y )Ax = 
j=l 1 2 2 j=l 2 2 

The indices i and j are as shown in Fig. 3. Upon rearrangement, one has 

with different notation 


V i,n+l/2 Ax^ j=l ^ U i+l/2,j u i-l/ 2 , 3 ^ 

where n is the number of fluid cells above the wall in the vertical 
direction. Note that (56) is not applied to fluid cells which are bounded 
on the top by the plate. The transverse velocity there is already forced 
to zero by the bound vorticity distributed along the plate. 

A final detail has yet to be mentioned. This pertains to the 
calculation of the diffusive flow of free vorticity from the top and bottom 
of the plate into the surrounding fluid. This involves the application of 
Eqns. (16), (17), and (18). 

The apparent slip velocity is found at each of 80 points on the 
plate, and yet there are only twenty fluid cells which are in contact with 
fluid on either the top or bottom of the plate. Therefore, the local slip 
velocity was integrated over each cell face, using Simpson's rule, in order 
to get an effective slip velocity for each fluid cell. This was not 
required on the wall since a local slip velocity was calculated at the 
center of each cell face in contact with the wall. 

Numerical Parameters. Most of the gross features of the 
computational grid are apparent from Fig. 2. There are 11,200 fluid cells 
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in all. There are 80 cells in the vertical direction and 140 in the 
streamwise direction. 

The cell dimensions are first determined on the plate. There are 
20 cells which span the unit plate length. First, 20 points are located 
over the unit interval according to the scheme in the 20-point Gauss 
quadrature formula. Control faces of fluid cells are then located midway 
between these points. This produces 20 cells with a size distribution 
which is symmetrical about x = 0, and which allows the clustering of cells 
near the plate leading and trailing edges. The same distribution is used 
in the fluid region above and below the plate, as well as ahead of and 
behind the leading and trailing edges. Outside the range — 1 < x £ +1 , 
the spacing is uniform and equal to Ax = 0.12 until x = +4.0. For 
x > 4.0, the horizontal spacing is 0.16. 

The unit length between the plate and wall has the same 
distribution of cells as that over the plate. Above the plate, cell 
dimensions expand in a fashion symmetrical to those below the plate, but 
for five cells only (rather than the usual ten). Beyond this point, the 
spacing is uniform and equal to Ay = .0517. 

The maximum time step was selected to be At = .01. Recall that 
this means that a fluid particle moving with the free-stream velocity will 
travel one plate length in 100 time increments. Provisions were made to 
reduce At, if necessary, in order to render the calculations stable. The 
reduction was always 10% of the current At. The At was never smaller 
than .0081. 
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Model for the Large-Scale Eddies 


A total of eleven different eddies was specified to enter the 
upstream flow boundary over the total time interval t = 10, For t > 10, 
the eddy pattern repeats itself. 

A schematic of the eddy pattern is shown in Fig. 4. The initial 
vorticity strengths are shown for each eddy. One can envision that this 
pattern moves to the right undistorted with the free-stream velocityjU^. 
The computational domain is just to the right. Therefore, this pattern can 
be supposed to exist for a distance of ten plate lengths upstream of the 
inflow boundary, and it moves as an ensemble until the outline of an eddy 
just touches the inflow boundary. After this point in time, the vorticity 
of the eddy, plus the background flow, convects into the domain with the 
local fluid velocity. 

Since the grid is rectangular and the eddy is taken to be circular, 
its shape could only be approximated. An example of this is shown for the 
third eddy in Fig. 5. The eddy spans vertical layers j = 7 to j = 17. 
The diameter is 0.7245, and the center is at y = .576 and l = 2.30. 

As seen in Fig. 5, horizontal grid lines intersect the circular 
outline of the eddy (dashed circle), thereby defining the sector of a 
circle. The difference between the areas of two sectors formed by adjacent 
grid lines corresponds to the area of a horizontal strip intersecting with 
the circle. A rectangle of equal area can then be formed by dividing this 
latter area by the height of a strip. In this way, a stair-step outline of 
the eddy is produced. 
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Fig» 4 Schematic Diagram of Eddy Pattern Introduced at 
Upstream Boundary 






Fig. 5 Representation of Circular Eddy by 
Rectangular Grid 
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The initial vorticity of each eddy is preassigned. The vorticity 
is assumed to be uniformly distributed over its area. The vorticity 
strength is chosen so that the velocity induced at the wall by a single 
eddy is no more than 2% of the free stream velocity. The nondimensional 
eddy strength varies from *128 to .362. To obtain the maximum perturbation 
velocity due to a single eddy, one divides the eddy vorticity by 4. 

The perturbation vorticity of the eddy is added to the background 
vorticity at the upstream boundary. This then defines the vorticity which 
enters an upstream control volume at any elevation, j. Vorticity enters 
according to a certain time schedule. Because the eddies are assumed to 
move toward the upstream boundary at uniform velocity U , and this 
velocity is used to nondimensionalize the time, their nondimensional 
locations ahead of the upstream boundary are equivalent to delay times. 
For example, suppose j = 10, and we are concerned with vorticity entering 
the grid from eddy #3 (shown in Fig. 5). For j « 10, the length of the 
equivalent rectangular area is 0.6865. The center of the eddy is at 
^ “ 2.30. Therefore, vorticity from this eddy is added to the background 
vorticity at j - 10 for ^ - t - t , where t = 2.3 - 0 .6865/2 and 
t 2 * 2.3 + 0.6865/2. 

Computational Procedure 

The calculations start from an initial state of rest. We envision 
that the flow is started impulsively with the background velocity field. 
The coherent eddies do not enter the upstream boundary until after a finite 
delay time. 
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The flow is initiated at t = 0 . Immediately, there is an apparent 
velocity slip along the plate given by < u slip > + p i ate = (D slip ) plate = 
U b (y=l). These apparent slip velocities form surface boundary conditions 
along the plate for the vorticity transport equation. As discussed in 
connection with Eqn. (32), the vorticity in the fluid cells adjacent to the 
plate is now non-zero and is given by +17^ (y= l) depending on whether the 
cell is above (minus sign) or below (plus sign) the plate. The vorticity 
throughout the entire field is calculated. A small lateral diffusion of 
background vorticity takes place. 

Note that the vorticity obtained in this first step includes the 
background plus perturbation vorticity. This will be the case for all 
calculations of the vorticity field. Also, the vorticity is convected by 
the total velocity (background plus perturbation). 

To obtain the perturbation velocities, however, we need only the 
perturbation vorticity, co *♦ This is obtained next by subtracting ft b (y) 
from the vorticity field. Now the boundary value for the perturbation 
stream function can be obtained from Eqns. (20) and (21), following which 
tjj 1 can be calculated. One contribution to the x-component of perturbation 
velocity is next obtained from 3^ f /3y everywhere in the field. 

It remains to calculate the per turbation velocities due to the 
bound vorticity, y, and the apparent slip velocities along the wall and 
plate surfaces. First, a) 1 is used to calculate the right-hand side of 
Eqn. (14), after whichy (x^) is obtained numerically from this integral 
equation. Finally, the contribution to the x-component of velocity due to 
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Y is obtained from the second integral in (12). Recall that the first 
integral has been replaced by 9^79 y. Now that both contributions to the 
perturbation velocity in the x-direction are known, these are added and the 
y-component of perturbation velocities is calculated from the continuity 
equation. 

The apparent slip velocities at solid surfaces are obtained from 
Eqns. (15), (16), and (17). These now form the boundary conditions for the 
vorticity transport equation. 

The computational cycle is now repeated. With each integration of 

the vorticity transport equation, the conditions at the upstream boundary 

are checked to see if eddy vorticity is ready to enter the grid. If not, 

then f^(y) is convected in. Otherwise, the eddy vorticity is added to 

ft , (y) and the total is convected in. 
b 


IV. RESULTS ABD DISCUSSION 


All computations were performed on the CDC CYBER 203 at the Langley 
Research Center. The computer jobs were submitted by remote batch from the 
Aerospace and Mechanical Engineering Department of the University of 
Arizona. The total central processing time for the results shown was 2.64 
hours . 

The primary results are presented in the form of vorticity contour 
plots. In addition, a sense of the flow development is provided by figures 
obtained using the marker and cell technique developed by Harlow and co- 
workers at the Los Alamos Scientific Laboratories (see Ref. [8]). 

The plots are grouped at the end of this section and are given 
after each 160 time steps, which corresponds to a dimensionless time 
interval of approximately 1.30. The first plot in each series shows 
contours of constant perturbation vorticity. This is the vorticity 
obtained after the contribution from the Blasius velocity profile has been 
substracted out. Contours are shown for the following values of 
perturbation vorticity, -.05, -.10, -.20, -.25, -.30, -.50, -1.00, + .50, 
+1.00. Following each vorticity plot is the flow pattern for the same time 
obtained with the marker and cell flow-visualization technique. It is 
emphasized that markers and cells were used only to help visualize the flow 
development after the complete flow field was obtained from the governing 
equations given in Chapter III. 
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The flow is started impulsively from an initial state of rest. One 
can see in Figs. 1 and 2 that the wake formed by the plate is developing at 
the same time that the vortical flow introduced at the upstream boundary 
( x = -i!u0 ) i s proceeding toward the leading edge of the plate. In Fig. 3, 
the eddies are just reaching the leading edge of the plate. Four eddies 
seem to be visible. Secondary vorticity is also being generated on the 
wall. The region close to the wall ahead of the plate contains positive 
perturbation vorticity, whereas that region below the plate contains 
negative perturbation vorticity. Although not shown on the figures, the 
wall region behind the plate contains slightly positive perturbation 
vorticity for a large extent behind the trailing edge. 

The regions of positive perturbation vorticity near the wall act to 
reduce the wall skin friction, since it combines with the negative 
vorticity of the background flow to produce a smaller velocity gradient. 
This indicates that there is an unfavorable (positive) pressure gradient at 
the wall ahead of and behind the plate. Directly below the plate on the 
wall there is a favorable (negative) pressure gradient. Clearly, the flow 
is accelerated as it flows between the wall and plate, but it is retarded 
ahead of and behind the plate. The local acceleration is caused by the 
thickening boundary layer on the plate which acts to restrict the flow area 
between it and the wall. 

The developing wake and flow perturbation ahead of the plate is 
also visible in Fig. 8A. Initially, all the flow markers were arranged 
uniformly over the flow field. 
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At the upstream boundary, new markers are introduced at a constant 
rate. This simulates hydrogen bubbles released uniformly from a wire 
placed across the flow. Because the fluid in the boundary layer (y < 2.0) 
moves more slowly than the outer flow, the markers tend to bunch up near 
the wall. One can also see that the particles appear to follow a wavy path 
ahead of the plate. This is because they follow the rotary motion due to 
the eddies. 

At an intermediate time of t = 8.14 shown in Figs. 11 and 11A, the 
eddies have merged with the boundary layers on the plate. Note that the 
perturbation vorticity of the eddies and wake has traveled approximately 
eight plate lengths, which is in accord with the elapsed time. The 
vorticity does not appear to penetrate the region directly below the plate, 
although it is prevalent below the elevation of the plate in the upstream 
region. This is probably due to the fact that the positive vorticity of 
the boundary layer on the underneath side of the plate is much stronger 
than the eddy vorticity. Thus when they merge, they combine to give a net 
positive vorticity. The vorticities of the eddies and boundary layer are 
of the same negative sign above the plate. 

Note also in Figs. 11 and 11A that there is a fine-scale eddy 
pattern in the far-wake region behind the plate. Furthermore, the negative 
vorticity adjacent to the wall directly below the plate has disappeared 
(compare with Fig. 9). This indicates that the enhanced skin friction 
there at earlier times has been reduced, as has the negative pressure 
gradient. Probably more of the flow is passing over the plate rather than 
between it and the wall. 
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For t > 8.14, only vorticity plots are given. The overall pattern 
follows that at earlier times. It is interesting that at t * 10.73 and 
12.03 (Figs. 13 and 14), the vortical pattern characteristic of "cat eyes" 
is swept over the leading edge of the plate with little distortion. 
However, at the last time shown of t = 14.61, the pattern is completely 
obscured. Clearly, the plate is acting to break up the eddies at the same 
time that it is preventing them from penetrating the region close to the 
wall. 

It appears from these computed results that the basic mechanism 
described by Corke, Nagib, and Guezennec [1] is correct. The plate 
prevents the large scale eddy from penetrating the wall region below and 
behind the plate. The strong wake produced by the plate blocks vertical 
excursions of high speed potential fluid into the boundary layer. Thus the 
vortical wake behind the plate is very straight and parallel to the wall. 
The eddies so prevalent ahead of the plate do not persist behind the plate. 
There does not appear to be a mechanism for their reforming into coherent 
structures. They seem simply to merge with the plate boundary layers. 

Before concluding this section, it is of interest to examine the 
time development of the drag on the system. This was computed at each time 
step from the momentum integral relationship. The contribution to the drag 
from the pressure field was neglected. The vertical flow boundary on the 
upstream side of the control volume used in the momentum balance was 
located at x — —3.4, That on the downstream side was located at x = +11,2. 
The drag force per unit of span is divided by the horizontal length of the 
control volume (14.6) to get the drag coefficient plotted in Fig. 17. 
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The drag coefficient is of limited quantitative value since it is 
not known what would be the drag if the plate were absent. However, some 
of the features do appear to be correlated with the actions of the eddies. 

There are a number of sharp peaks in the drag curve. These occur 
at approximately t = 1.3, 2.5, 4.2, etc. It appears that these peaks are 
associated with the arrival of an eddy core at the upstream boundary of the 
momentum control volume. Since the eddy pattern repeats itself for 
t > 10, the peaks at t = 11.3 and 12.7 appear to be repeats of those 
occurring at t = 1.3 and 2.5, although they occur at a higher level. 
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Fig, 6 Contours of Constant Perturbation Vorticity. 
t = 1.44 
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Fig. 8 Contours of Constant Perturbation Vorticity. 
t = 4.25 
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Fig. 15 Contours of Constant Perturbation Vorticity. 
t = 13.33 
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Fig« 16 Contours of Constant Perturbation Vorticity. 
t = 14.61 
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Fig. 17 Development of the Overall Drag Coefficient with Time 
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V. SUMMARY AMD CONCLUDING RKMARTTS 


Numerical solutions of the unsteady Navier-Stokes equations have 
been used to simulate large eddy breakup by a single flat plate imbedded in 
a wall boundary layer. The flow field has been treated as two-dimensional 
and incompressible. 

The plate length has been taken to be one— half the boundary— layer 
thickness. The Reynolds number based on plate length was held fixed at 
1000, giving a Reynolds number based on boundary-layer thickness of 2000. 

The numerical simulation is essentially that of a laminar boundary 
layer undergoing transition to a turbulent boundary layer. Vorticity 
perturbations in the form of coherent eddies are swept into the 
computational domain. All of the eddies interact with one another and the 
background flow in a fully nonlinear fashion. 

The eddy structures are altered by the plate imbedded in the 
boundary layer. The strong vorticity layers produced by the plate merge 
with the vorticity of the eddies, and the resultant wake persists for a 
long distance behind the plate. The coherent structures so evident ahead 
of the plate are nonexistent behind the plate. The straight appearance of 
the wake region reveals that there are no strong transverse velocities 
behind the plate. This is not the case ahead of the plate. Therefore, the 
plate effectively straightens the slow behind it, and high speed fluid 
above the plate does not penetrate the region close to the wall. 
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It is concluded that the plate does suppress the lateral mixing of 
the fluid in the boundary layer. This is consistent with the explanation 
given by Corke, Nagib, and Guezennec [1] for one of the mechanisms 
contributing to the drag reduction revealed by their experimental 
measurements. Although their measurements were made in a more highly 
developed turbulent flow than that simulated here, the fundamental 
conclusion is the same. 

The present numerical simulation leaves unanswered the question of 
whether or not the plate actually reduces the total drag. To partially 
answer that question, one would need to perform a companion calculation 
with the same flow parameters, but with the plate absent. A comparison of 
the resulting drag curves would then show if the presence of the plate 
actually reduces drag. 

Whereas the above-mentioned comparison would be of interest, it is 
beyond the scope of the current study. Also, the outcome would not be 
totally conclusive since many features of a turbulent flow are absent in 
the current simulation. The most important of these is the three- 
dimensional character of the flow. The introduction of a third dimension 
would allow the vorticity to change due to the stretching of vortex lines. 
That is, vorticity could increase locally. In the current calculations, 
vorticity only decreases due to the mechanism of diffusion. Local regions 
of enhanced vorticity cannot occur, and thus one mechanism for lateral 
mixing is completely absent. 
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APPENDIX A 


SOLUTION ALGORITHM FOR THE STREAM FUNCTION 


The algorithm is based on the Stabilized Error Vector Propagation 
method, as described briefly in Section III and Ref. [4], The method is 
based on the observation that if the stream function is known on two 
adjacent rows, then it can be advanced to the next row. In this way, the 
solution can eventually be advanced to all the rows in the computational 
domain. 


The governing equations for ip are given by Eqns. (34), (44), (45), 
(46), and (47). We assume that 


- '■p. . + e. . 

y i,J y i,J i,j 


(A-l) 


/\ 


where ip. . is the exact value, and E, . is an error term. Eqn. (A-l) is 
substituted into (37), (47), (48), and (49). The result is a set of 
algebraic equations for £. . For purposes of illustration, we assume a 
uniform grid with Ax^ = Ay \ . Then one has 


e. . , = 4e. .-e.,, . - e. . . - e. ... 
i,j-l i,j i+l , 3 i-l»J i.J+l 


(A-2) 


for i f 1, 141 and j f 80. For the nodes adjacent to boundaries, we 
have to make the following replacements, 
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3 


(A-3) 


_ 4e. . - e.., . 

■ .. , — i±ia 


P _ 4e. . - e. . . 

* 1 + 1.3 ' 


E i,j+1 - £ i,j 


(A-4) 


(A-5) 


Notice that the exact boundary conditions are enforced on the sides and top 
of the grid. 

The main task is to trace the propagation of errors introduced at 
the top of the grid to the bottom of the grid, where a final boundary 
condition can be enforced. For example, consider the calculation 
illustrated in Table 1. 


Table 1. Illustration of Error Propagation. 


es 


2 

3 

4 

5 

6 

7 

80 

0 

0 

1 

0 

0 

0 

0 

79 

0 

-1 

3 

-1 

0 

0 

0 

78 

2/3 

-7 

13 

-7 

1 

0 

0 

77 

52/9 

-122/3 

63 

-41 

11 

-1 

0 

76 

1130/27 

-2020/ 9 

962/3 

-231 

85 

-15 

1 


A unit error introduced at i = 3, j — 80 propagates through the grid as 
shown. At j =76, the error is already quite large. 
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Note that the error at any j is directly proportional to the error 

introduced at i = 3, j = 80. Let us first focus our attention on the error 

coefficient at i s 1, j = 76. In this case let us call it a ^3 = 1130/27. 

The first subscript denotes the column in which the error coefficient 

occurs, and the second denotes the column in which the unit error is 

introduced. If an error e is introduced at i = 3 and j — 80, then the 

3 

error at i - 1, j = 76 is e 3 * 

Now observe that if the unit error had been introduced at i = 2, 
rather than i = 3, a different number would have occurred in the first 
column (i = 1) at j = 76. Call this the error coefficient, a 12 * Likewise, 
we could generate different coefficients a^ in the first column at j -76 
for each unit error introduced at column i. The total error at node i = 1, 
j = 76 is thus 


'1,76 “ a ll e 'l + a 12 e 2 + a 13 e 3 + * 


(A-6) 


where only has actually been computed in this example. 

To obtain the error at j = 76 and i = 2, we need to generate 
coefficients a 2V a 22 , a 23 , etc. In this example, we have obtained only 
a 23 = -2020/9, but one would have in analogy to (A-6) 


E 2,76" a 21 e l + a 22 e 2 + a 23 e 3 + ’*• 
and in general, one can write 

E n, 76 = a nl e l + a n2 e 2 + a n3 e 3 + "* 


(A- 7) 


(A- 8) 
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Thus we see that the a n ^'s form the error coefficients of an error 
propagation matrix. The complete error vector at Row 76 is just 

tE 76 ] = [A] [e] (A- 9) 

If we had wanted the error at Row 77, then we would have obtained a 
different error propagation matrix, [A], 

The point to (A- 9) is that we can relate the error introduced at 

one row to the error at any other row. For example, we guess the value of 

the stream function along the top row j = 80. For this study, it is 

assigned the arbitrary value of zero. The stream function is calculated by 

means of Eqns. (37), (47), (48). (49) until we reach the row along the 

wall, for which ^ is given fay (50). Now this row of values should all 

be zero, according to the boundary condition required at the wall. 

However, it won't be zero because an incorrect guess was made for ^ ™ 

1 , 80 . 

In fact, the values of f comprise the final error vector at the wall. 

1 , w 

Call this [E w ] . We then have by (A-9) 

[E w l = [A] [e] (A- 10) 

] and [e] are vectors with 140 elements. [A] is a square matrix. Since 

IV and tA] are known, we solved for [e] to get 

[e] = [A]~l [EJ (A-U) 
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Thus we have the error vector along the first row of nodes. This can be 

used to correct the values of ^ . QA , which were initially guessed to be 

1, oU 

zero. Once the exact values of . OA are known, the true solution for the 

1 , oU 

stream function can be recalculated throughout the whole domain. 

Thus there are two passes through the calculations, once the error 
coefficient matrix, [A], is generated. The first pass produces the error 
vector at the wall, and the second pass produces the stream function. 
Therefore it is a direct solution method. The advantages of this method 
are twofold. It is fast because it is noniterative. Also, the maximum- 
error can be found. This error is the deviation from zero of the final 
stream function value along the wall. Typically, the error is close to the 
computational precision of the computer. 

There is a major disadvantage, however. The error matrix [A] can 
become so ill-conditioned that [A]”l is inaccurate, and thus the error 
calculated from (A-ll) is inaccurate. This happens when there are many 
rows in the grid. To handle such cases, the grid has to be broken down 
into blocks. This allows the error to be stabilized. 

The algorithm used for the Stabilized Error Vector Propagation 
(SEVP) method is illustrated here by the use of three blocks, as shown in 
Fig. A-l. It can be generalized to any number of blocks. Note that 
adjacent blocks overlap by two rows. The values for the stream function 
are guessed on dashed lines. Boundary values for the stream function are 
given on solid lines. Suppose each block is 8 rows long. We start at the 
top of the first block. 
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i=l40 


_ [«>]_ 

[ e 2] * [A||][e|] 

[e 3 ]» [A|2][e,] 

[*4]; i A _2lK e j]_ 
[ e 5 ] 1 [ A 22l( e 3] 

C e 6l* [ A 32K e 5l 


BLOCK 2 


tration of the Use of Blocks in the 
lized Error Vector Propagation Method 





Error propagation matrices [A^], [A 12 ], [A 2] .], ^22 ^ ’ and t A 32^ 
must first be generated. This can be done once and the results stored 
for later use. We first specify 140 unit vectors on j = 21, j = 14, and 
j = 7. Recall that this is done by making an entry of unity at one of the 
columns i; all other entries being zero. We then calculate the errors by 
Eqns. (A-2) through (A-5) at each of the j rows. This requires "boundary 
values" for the errors to be specified at j = 15 and j = 8. Eqn. (A-5) is 
enforced for j = 21 only, since this is a true boundary for the grid. 

Calculations begin in the first block, and [A^ 2 ] and [A^ 1 
are obtained. Clearly, [e 2 ] - [Ajj] f e ^] and [e-j] “ [A^ 2 ][e^]. Then 
[e^ = [A^] - ! te^lj whereupon 

[e 2 ] = [A.^] (A- 12) 

We now introduce unit vectors along j = 14 and enforce values of [e 2 ] 
obtained from (A-12) along j = 15. We thus march to Rows j = 8 and j = 7 
and obtain [ and [A 22 ] • Once again 

[e^] = lA 21 l [e^l = lA^] [A 22 ]“l le 5 ] = [ B 2 ] [ e 5 ] (A-13) 

Unit vectors are introduced along j - 7, and corresponding boundary values 
[e^] are calculated from (A-13) and enforced along j = 8. The solution for 
the error vectors is again extended into the third block to j = 0. This 
gives IA 32 ]. 

Error propagation matrices [A^l, [A^ 1 > [A^l, ^‘22^ and ^A^], 
and corresponding matrices [B^] and [B 9 ] are now all known and stored once 
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and for all* It remains to find an approximate solution for the stream 
function. This requires two forward sweeps through each of the blocks. A 
final backward (correction) sweep through the entire grid produces the 
exact solution. 

We start by guessing the stream function to be zero along Row 
j = 21 and use the boundary condition at j = 22. We march the solution 
forward to j ■ 15 and 14. Along j = 14, we must enforce a pseudo- 
boundary condition, which we choose to be ip = 0. The departure from zero 
of the calculated value of ^ on j = 14 gives £e J. 

Next, we obtain [e^] = [©gl* With this we correct the 

solution for ^ along j = 21 , and we generate stream functions over the 
first block and into the second block until we reach j = 7. We again 
enforce the pseudo-boundary condition of ^ = 0 here, and obtain the error 
vector [e 5 l. From this we find [ e 3 ] = [ A 22 ] “* 1 [ e 5 3 and t e 2 ] » [B^ [ e 3 ] • 
With these boundary values, we now correct the stream function along j = 15 
and 14 in the second block and proceed to march ^ through the second and 
third blocks. The process is repeated until we have swept all the blocks. 

For the last block, we use the exact boundary condition of ^ = 0 
on j = 0 to find [e 6 ]. Then [ e 5 ] = [ A 32 ] “1 [ e g ] and [ e 4 ] = [ B 2 ] [ e 5 ] . The 

solution for ip is corrected on j — 7 and 8 and marched through the last 
block to j = 0 . 

At this stage we have the exact solution for ^ over the last 
block. It remains to extend the solution, along with new error vectors, 
over the second and first blocks. For this we obtain a new error vector 
[e^] by subtracting the exact solution for ^ at j = 7 from the approximate 
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solution generated at j = 7 during the forward sweep of Blocks 2 and 3. 
Note that this approximate solution along j - 7 was stored with 
forethought. We now recalculate [e^] = tA 22 l~^ [e^] an< * = 

Then, the approximate solutions at j = 14 and 15 (again stored with 
forethought) are corrected and the new exact solution for ’P is marched 
from j = 15 to j = 7. A new error vector [e^] is obtained by subtracting 
the exact solution on j = 14 from the approximate solution stored at j = 14 
from the forward sweep through Blocks 1 and 2. Finally, [e-^] = [Aj^ ]"^ 
[e^], and the approximate solution on j = 21 is corrected. The boundary 
condition along the top of the grid is applied to find rj» along j = 22. 
The final sweep from j = 22 to j = 14 completes the generation of the exact 
solution. 

The solution generated by this direct method is discontinuous at 
block boundaries (j = 14 and 7), and it may not satisfy the exact boundary 
condition of $ = 0 along j = 0. The amount of discontinuity and the 
departure from 4* = 0 on j = 0 gives a measure of the maximum error. 
Typically, the maximum error can be controlled to the round-off error of 
the calculations. Iterations can be used to reduce errors further, but 
they were not used in this present study. 

In this work, fourteen blocks were used. The first ten blocks had 
eight rows each, and the remainder had seven rows. The maximum error on 
j = 0 was 1.1 x 10“li at i = 36. Most of the errors were on the order of 
10“ 15 along j = 0. Errors at the block boundaries were not checked. 
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APPENDIX B 


INTERPOLATION METHOD FOR STREAM FUNCTION 

The interpolation procedure involves a nine— point formula in 
x and y. To proceed, we assume a biquadratic formula of the following 
type: 

^ — a + bx + cy + dxy + ex2 + fy2 + gx2y + hxy2 + ix2y2 (B— 1) 

Here, x and y denote local coordinates, with origin at Point 1, as shown in 
Fig. B-l. 

In principle, one can find the nine unknown constants a through i 
from known values of \p at nine points in the field. However, the solution 
of the nine simultaneous equations would be tedious, and we do want a 
closed form solution. Instead, we utilize a procedure used in finite- 
element analysis. 

We note that the x— and y— spacing is variable and construct 
polynomials in x as follows: 


(x-h^ ) (x-h^-h^) 


N l(x) 

1 ' ' 1 2 
h 1 (h ;] +h 2 ) 

- x(x-h.-h 0 ) 

(B-2) 

N2< x ) 

_ 12 

h l h 2 

x(x-h^) 

(B-3) 

N 3 (x) 

- (h 1+ h 2 )h 2 

(B-4) 
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These expressions 

have 

the following 

properties : 

N^O) 

- 1 

* 

N 1 (h 1 ) = 

N l (h l +h 2 } = ° 

N 2 Ch L ) 

= 1 

f 

n 2 (o) = 

N 2 (h 1 +h 2 ) = 0 

N 3 Ch l +h 2 ) 

= 1 

> 

n 3 C0) = 

N 3 (h l +h 2 ) = 0 


Next, we form identical polynomials in y, but with replaced by 
k^, etc. The final expression for ip becomes 

ip = N^Cx) + to 2 N 2 (x) + N^Cx)] N (y) 

+ N 1 (x) + w 5 N 2 (x) + N 3 Cx)] N 2 (y) 

+ ^7 N i^ x ^ + u 8 N 2^ + “9 N 3^ x ^ N 3 (y) (B-5) 

Next, Eqn. (B-5) is expanded, and terms are collected. After extensive 

algebra, one can obtain the coefficients a through i by inspection. 
The results are: 


a = ^ (B-6) 


-(2h 1 +h 2 ) 

h l+h 2 

h l 


h 1 (h 1 + h 2 ) 

“1 + h x h 2 “2 

h 2 (h 1 + h 2 ) W 3 

(B-7) 

-( 2k x +k 2 ) 

k l +k 2 

k l 

(B-8) 

k^k^ k 2 ) 

1 k 1 k 2 w 4 

k 2 (.k 1 + k 2 ) w 7 
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d = 


Uh-j+hp (2k 1 +k 2 ) ( h ! + V ^ 2k 1 +k 2 ) 

h l k l Ch l +h 2 ) (k l +k 2 } ° J] - " h l h 2 k l (k l +k 2 } ^ 


h 1 C2k 1 +k 2 ) 


k l h 2^ h l +h 2^ ( k l +k 2^ 


CO- 


(2h 1 +h 2 ) (k 1 +k 2 ) 


w, 


h l k l k 2^ h l +h 2^ 4 


(h 1 +h 2 ) (k 1 +k 2 ) 
\^h~ 2 k^ W 5 


h 1 (k 1 +k 2 ) 


CO, 


k 2 k l k 2^ h l +h2 ^ 6 


k 1 (2h 1 +h 2 ) ^ 

h^k 2 (h^+h 2 ) (kj+k 2 ) 


7 - 


(h 1 +h 2 )k 1 "1"1 » 

h ] _h 2 (k 1 +k 2 ) 0 ^ h 2 k 2 (h 1 +h 2 )(k 1 +k 2 ) 


h l k l 


e = 


hi'CVV Wl ” h l h 2 ^ + VW ^ 


f = 


k 1 (k 1 +k 2 ) W 1 " k 1 k 2 “4 + k 2 (k 1+ k 2 ) “7 


g = 


- ( 2k^+k 2 ) 


h 1 k 1 (h 1 +h 2 )(k 1 +k 2 ) “l + h^^Ck^kj,) “2 " k 1 h 2 Ch 1 +h 2 )(k 1 +k 2 ) 


2k 1+ k 2 


co„ - 


2k x +k 2 


k 1 +k 2 


k l+ k 2 k 1 +k 2 

co, - ; — co r + 


co. 


+ h ] _k 1 k 2 (h 1 + h 2 ) 4 h 1 h 2 k ] _k 2 h 2 k 1 k 2 (h 1 +h 2 ) 6 


h ] _k 2 (h 1 +h 2 )(k 1 +k 2 ) “7 + h 1 h 2 k 2 (k 1 +k 2 ) “8 h 2 k 2 (h 1 +h 2 (k 1 +k 2 ) “9 


h = 


-(2h 1 +h 2 ) + “i‘“2 - “1 c o 

h 1 k 1 (h 1 +h 2 ) ( k i +k 2 )_ h 1 h 2 k 1 (k 1 + k 2 ) “2 k i h 2 (h l +h 2 ) (k l +k 2 } 3 


h^+h 2 


2h_+h 0 h-+h- 

1 2 G), 1 2 0)- 

+ h 1 k 1 k 2 (h l +h 2 ) 4 " h ! h 2 k i k 2 hjc,k 0 (h,+h,) 


<0, 


2 12 v 1 2 J 


2h x +h 2 


h l +h 2 


h l k 2 Ch l +h 2 )(k l +k 2 )W7 + h l h 2 k 2 (k l +k 2 ) ^ h 2 k 2 Ch 1 +h 2 )(k 1 +k 2 ) “9 


(B-9) 

(B-10) 

(B— 11) 
J 3 

CB-12) 

CB-13) 
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Oj, 


CO. 


CO. 


X = 


h l k l^ h l +h 2^ k l +lc 2^ h^t^k^k^+k^ k^h 2 Ch^+h 2 > (k^+k 2 > 


(O, 


co. 


co. 


h l k l k 2 ( ' h l +h 2' ) 


CO- 


k l h 2 k l k 2 


co 


8 


h 2 k l k 2^ h l +h 2^ 


co. 


h l k 2 (h l +h 2 )(k l +k 2 ) h,h 0 k 0 Ck,+k 0 ) h n k 0 (h 1 +h 0 )(k 1 +k 0 ) 


(B-14) 


12 2 v 1 2 


2 2 v l 2 ' x 1 2 


The spacings h_j, , etc. are obtained from the nodal arrangement 
shown in Fig. 3; that is, h-^ = .5 (Ax^_^ + Ax^) etc. 
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